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1. Introduction 

While QCD is accepted as the correct theory of the strong interactions, it is frustrating that we 
still are unable to derive from it the basic properties of matter or map out a phase diagram that is 
more than a cartoon. The problem is that the only tool we have for analyzing QCD nonperturba- 
tively outside of special kinematic regions is the lattice, and lattice QCD faces daunting barriers 
to investigating systems with significant baryon number. In the grand canonical formulation, the 
problem encountered is that the fermion determinant has a rapidly varying phase leading to strong 
cancelations in the path integral and a suppression of the partition function that is exponential in 
the volume. This is often called the "fermion sign problem", although we will argue it is not spe- 
cific to fermions. In a canonical formulation, where one computes correlation functions of fixed 
baryon number, one encounters a couple of problems. The first is that the measurements may be- 
come increasingly noisy as the baryon number is increased, which we will call a "noise problem". 
The second is that the value obtained for the correlator may be relatively quiet, yet drift with time, 
giving no evidence for a plateau in the effective mass from which one can infer a ground state 
energy. The effective mass computed from a correlator in a system with discrete nonzero energy 
levels must exhibit a plateau at late time; so if a plateau is absent it follows that the path integral 
is not being estimated correctly, with little overlap between the field configurations being sampled 
and the ones that actually dominate in a correct estimation of the correlator. We will refer to this 
second problem as an "overlap problem". In an actual simulation both problems may be found 
simultaneously: a drifting and noisy value for the correlator. 

After arguing that the problems encountered in the canonical and grand canonical approaches 
arise from the same physics, we focus on the canonical method and the statistical distribution of 
correlators. We give examples where the overlap problem is a manifestation of an underlying 
statistical distribution for the correlator which is very heavy-tailed, and not too far off a log-normal 
distribution. We show that the appearance of nearly log-normal distributions is indicative of a mean 
field theory expansion, such as that used to study electron properties in disordered media. We then 
develop a technique to effectively overcome the overlap problem by defining mass plots using an 
estimator which is not minimum bias, but where the log-normal result is the first term in a well- 
defined expansion. This method has been used with great success in studies of trapped "unitary 
fermions" — a strongly coupled nonrelativistic conformal theory relevant to trapped atoms at a 
Feshbach resonance — allowing a measurement of the ground state energy of up to 70 fermions in 
a harmonic trap. We also show QCD data from the NPLQCD collaboration that indicates similar 
techniques might aid in the study of QCD; we expect this to be especially so at larger baryon 
number. This work has been published in Ref. [1]. 

As a possible direction for future research, our approach suggests a renormalization group ap- 
proach to the statistical distributions of lattice observables, and a way to analyze such distributions 
in a manner analogous to effective field theory, where the n^^ cumulants of the log of the observable 
play the role of higher dimension operators in the effective Lagrangian. 

2. Noise and the physical spectrum 

We first examine the connections between the noise problem encountered in canonical sim- 
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Figure 1: Effective mass plot for nine quarks in the triton channel from NPLQCD [7]; annotation by D.K. 



ulations and the phase problem due to the fermion determinant in grand canonical computations, 
developing a heuristic picture for why a sign/noise/overlap problem arises in the first place. 

Grand canonical simulations of lattice QCD at nonzero quark chemical potential involve a path 
integral measure that includes a fermion determinant which is neither real nor positive, leading to 
an exponentially hard computation [2]. It was noticed long ago [3] that the determinant phase 
problem sets in for quark chemical potential li > mji/2, even at zero temperature, which seems 
peculiar since nothing physical happens until /i > itin/^, being the nucleon mass. A beautiful 
explanation for this precocious annoyance was given by Splittorff and Verbaarschott [5] (see also 
[6] for earlier related ideas): they pointed out that in two flavor QCD with degenerate quarks, the 
quark determinant phase is exactly what sets apart a simulation at finite isospin chemical potential 
(Mw = — = M) from a simulation at finite quark number {jiu ^ jid^ ji). In the former system, 
pion condensation occurs at /i = m;7r/2 leading to a drop in the free energy density, and therefore 
growth of the partition function Z which is exponential in the volume. Therefore a simulation at 
finite baryon number with ji >mji/2 must have the quark determinant phase oscillate wildly to 
cancel this exponential growth, since there is no pion condensate in this system. The sign problem 
in QCD is therefore characterized by A/i = ^ = (Mn/3 — mji/2), the difference between the values 
for ji where the ground states rearrange themselves in the isospin and baryon number systems. 

In a canonical simulation of the A-nucleon state one looks instead at the 3A quark correlator 
over Euclidian time T, Ca(t,(7), which is a function of the link fields U for the gluons. Eor short 
times, the correlator gets contributions from the many excited states with 3A quarks, but at later 
time it should be dominated by the ground state for A nucleons. One typically creates an effective 
mass plot of 

meff(T) = --^ln(CA(T + AT,C/))/(CA(T,C/)) (2.1) 

At 

where the average is over an ensemble of link fields and At = 0(1), chosen to optimize the result. 
When the correlator is dominated by the ground state, the effective mass should be a constant, the 
ground state energy. Eig. 1 shows an effective mass plot for the case for A = 3 provided by the 
NPLQCD collaboration clearly exhibiting three types of behavior: strong T-dependence at short 
times when excited states contribute, a plateau dominated by the mass of the triton, and degradation 
into noise at later time. 
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The growth of noise in the correlator makes it difficult to extract a signal without a lot of 
statistics, and the problem is expected to get worse with increasing baryon number. It would have 
been strange if properties of matter at finite baryon density, which are extremely difficult to compute 
in the grand canonical ensemble, were easy to compute in a canonical simulation — and indeed this 
noise problem exhibited by the correlator at late time is presumably an avatar of the sign problem. 
This sign/noise problem does not arise simply because of Fermi statistics; for example, constructing 
correlators Ca as 3A x 3A Slater determinants of quark propagators accounts for Fermi statistics 
and leads to a computational cost from the determinant only scaling as (3A)^, not the exponential 
difficulty seen in the grand canonical computation. Instead, the noise is closely related to the 
physical spectrum, as has been quantified by Lepage [4]. For example, in QCD the expectation of 
a 3A quark correlator for a nucleus of atomic number A and mass Ma is (Ca) ^ e~^^^, while the 
variance in the sample mean Ca can be estimated as 

= -1 ((CIQ) - (C1)(Q)) ~ -1.-3^-^ (2.2) 

for sample size Since Ca corresponds to 3A quark propagators and C| to 3A anti-quark prop- 
agators, the variance is dominated by the state with 3A pions. Thus the signal to noise ratio scales 
as ^ exp(— 3A^t), where ^ = {Mn/3 — mji/2) is the same parameter we saw characterizing 
the sign problem in the grand canonical case. This reasoning is rather simplistic, as the overlap 
between the operators and the nucleon or pion states — which will typically contain volume fac- 
tors - has not been included. However, Fig. 2 shows evidence from QCD simulations by NPLQCD 
that the Lepage argument is qualitatively correct. The noise and sign problems are therefore pre- 
sumably closely related and determined by the physical spectrum of the theory and should not be 
thought of as solely a "fermion sign problem"; similar issues can also arise in interacting boson 
theories. 

In either case the problem arises from the existence of multiparticle states for which the en- 
ergy/constituent is lower than for the states one wants to study: in QCD, quarks are in some sense 
lighter when they are in a pion than when they are in a nucleon, by the amount ^. This is a problem 
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Figure 2: A plot of — l/dn(signal/noise ratio) versus time for p (red) and pp (purple) correlators computed 
in lattice QCD; horizontal lines give the expected values for these quantities using the Lepage argument. 
Figure supplied by K. Orginos (NPLQCD), annotation by D.K. 



4 



Listening to Noise 



David B. Kaplan 



N=46 (Ni=N^=23), L=12, Ar=3 




Figure 3: The effective mass plot for 46 unitary fermions in a box of size L = 12. 



that can occur in bosonic systems as well, and so one should not expect the sign problem to simply 
be a fermion problem. Simple bosonic path integrals in the grand canonical approach, such as 
gauged theories with a 0^ interaction, manifestly have complex action in Euclidian spacetime and 
a sign problem. 

Applying the Lepage argument for Ca(t, (7) to higher moments one sees that all odd moments 
are proportional to exp(— M^t) while even moments are all dominated by pion states^ . Thus at very 
late time T, the probability distribution for Ca(t,(7) will be almost symmetric, implying a noisy 
sample with small mean. At earlier times, however, we will see that the distribution appears to be 
very asymmetric and heavy tailed in both QCD and other theories, and it is this regime we wish to 
analyze. 



3. Unitary fermions and a Mean Field Description 

Nonrelativistic spin half fermions with strong short-range interactions tuned to a conformal 
fixed point where the phase shift satisfies d{k) = k/2 for all k are called "unitary fermions". One 
can think of the inter-particle interaction as being described by a square well potential tuned to 
have a single bound state with zero binding energy, as one takes the width of the potential to 
zero. The field theory describing this system is conformal, and interesting to study both for its 
simplicity and universality, its challenges for many -body theory, and because it can be realized and 
studied experimentally using trapped atoms tuned to a Feshbach resonance [8]. It is also an ideal 
nonperturbative theory for studying fermion sign problems on the lattice, being simpler and faster 
to simulate than QCD. At its most basic, the lattice action involves the simplest discretization of 
the Euclidean Lagrangian [9] 

xi/\d^ - VV2M)i//- ^m^0^ + (3.1) 

where is a non-propagating auxiliary field with tuned to a critical value, and i// is a spin ^ 
fermion with mass M; a more sophisticated action tuned to reduce discretization errors was recently 
presented in [10, 11]. 

^Thanks to Martin Savage for first pointing this out to us. 



5 



Listening to Noise 



David B. Kaplan 



The ground state for N = (A/| + A/^) unitary fermions can be determined by computing the 
correlator for each background field, and then averaging over the ensemble of such fields. One 
finds that when N is large, the correlators rapidly become extremely noisy and drift away from a 
plateau; see for example Fig. 3 for an effective mass plot for N = 46 unpolarized fermions in a 
box with periodic boundary conditions of size L = 12. In order to better understand the problems 
apparent for t > 10, it is instructive to plot a histogram of the correlator before averaging over 
the field. This is shown in Fig. 4, which reveals a distribution for A/^-body correlators Ca^(t, 0) 
which is increasingly non-Gaussian at late t; in fact. Fig. 4 shows that it is In Cn which appears to 
be roughly normally distributed, so that Ca^ (t, 0) is approximately log-normal distributed with an 
increasingly large <7 and long tail at late time. 

A log-normal distribution for the correlator is difficult to reconcile with the Lepage argument. 
First of all, Lepage assumes that the central limit theorem applies and that the distribution of cor- 
relators is close to Gaussian and well described by the mean and variance, which is evidently not 
the case here. As we will see in a toy model below, the Central Limit Theorem may be irrelevant in 
such cases, requiring an exponentially large number of configurations to be applicable. Secondly, 
the size of the variance for this log-normal distribution cannot be related to the physical spectrum of 
the theory using the Lepage argument. The Lepage argument relates the variance in measurements 
of Cn to the energy of the lightest state coupling to (Ca^)^; there are no anti-particles in this theory, 
and and so {Cn)^ couples to the ground state of a system N particles each of two different species 
of spin 1/2 unitary fermions ^. However, such a theory in the continuum does not have a finite 
energy ground state, being unstable against condensing into a point-like clump. On the lattice, a 
ground state exists for such a system since there is a maximum density due to Fermi statistics and 
the cutoff; this ground state is a lattice artifact, however, and even then may have a vanishingly 
small probability to be created by (Cn)^, in which case the variance might be determined by some 
metastable or excited state. 

The appearance of a heavy-tailed distribution should not be surprising, however, since having 




log(c) 



Figure 4: Distribution histograms for c = C/v(t, 0) and ln(c) for A/^ = 4 unitary fermions at several times 
T, taken from Ref. [11]. Curves fitting ln(c) are Gaussian, implying that c is approximately log-normal 
distributed, with increasing with time. 



^If C/v is the product of A^^ antisymmetrized correlators times A^^ antisymmetrized correlators, then {C^)^ corre- 
sponds to a product of four groups of antisymmetrized correlators, and is equivalent to the propagator for two distin- 
guishable species of spin half fermions. 
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fermions wander through a random background is similar to the problem of electrons propa- 
gating in disordered media, where heavy-tailed distributions are ubiquitous in the vicinity of the 
Anderson localization transition. For example, it is found that for physical quantities such as the 
current relaxation time or normalized local density of states, the distribution function P(z) scales 
as exp(— Qln^z). A particularly simple way to derive these results is to use the optimal fluctuation 
method of Ref. [12], which is a mean field approach. We can adapt these methods to the current 
problem, defining the variable Y = IuCa^ (t, 0) and computing its probability distribution P{y) as 

P{y) oc I D^e-'<^ 5(F(t,0) -y) = J ^e'' (3.2) 



where S(j) = J d^x^cj)^ and S = S(j) —it(lnCN{T,(j)) —y). Using the PDS subtraction scheme [13] we 
have = MX /An, where the renormalization scale A is taken to be the physical momentum scale 
in the problem — in this case X = kf = {Sn^N/Vy^^, N/2 being the number of fermions with a 
single spin orientation. We proceed now to evaluate this integral using a mean field expansion; it is 
not evident that there is a small parameter to justify this expansion, but the leading order result is 
illuminating and fits the numerical data well. We expand about 0(x) = ^ = to, and use the fact 
that for large T the n^^ functional derivative of IuCa^ (t, 0) with respect to (x) equals the the 1-loop 
Feynman diagram with n insertions of Xf/^xi/ in the presence of a chemical potential ji^k]^/ (2M). 
The equations for 0o and are given by 

. m^0o .ym^0o 

y-lnZ + T£o(jV) 
</>o = (3.3) 

where Eo(N) = 3NEf/5 is the total energy of N free degenerate fermions (N/2of each spin), and 
Z is the overlap of the source and sink with the free fermion state. The leading term in the mean 



field expansion for P{y) can therefore be expressed as P(j) oc exp 



iy-yy 



with 



3; = InZ - TEoiN) , = ^^o(A^) T . (3.4) 

This describes a log-normal distribution for the A/^-fermion propagator CAf(T,0), with both mean 
and variance growing in magnitude with time in units of the energy of N free degenerate fermions. 
In Fig. 5 we plot the quantities — and as a function of N obtained from correlator 

distribution data for unitary fermions at late T, and find that the gross features of the results are 
compatible with the mean field estimates of unity and 40 /9n obtained from eq. (3.4). 



4. A toy model 

It would be useful to devise an algorithm to reliably estimate energies without having to ex- 
haustively sample the long tail of the correlator distribution, yet without making incorrect assump- 
tions about the exact functional form of that tail. An approach we suggest here is to exploit the 
general relationship between stochastic variables X and Y = InX: 

ln(X) = £^ (4.1) 

n=l 
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where is the n^^ cumulant of Y. This relation can be proved by noting that the generating 
function for the is ln07(r) where 0f(O = {e^^) = {X^) is the moment generating function for 
F, and evaluating at = 1, assumed to be within the radius of convergence. The motivation for 
investigating eq. (4.1) is that if the distribution P(X) were exactly log-normal, the above sum 
would end after the second term, as K'„>2 would all vanish; therefore by replacing the by sample 
cumulants and truncating the sum at n = n^ax, one might hope to have a reliable estimator for 
ln{X), provided that the systematic error from truncating eq. (4.1) and the statistical error from 
sampling Kn can be simultaneously minimized. 

Distributions with log-normal-like tails arise naturally in products of stochastic variables. The 
propagator Ca^ (t, 0) for unitary fermions can be expressed in a transfer matrix formalism as the 
product of T matrices — one per time hop — each of which is the direct product ofNVxV matrices 
of the form e~^/^{l + g(p)e~^/^, where is a constant matrix (the spatial kinetic operator), (p is a 
random diagonal matrix with 0(1) entries corresponding to stochastic fields living on the time 
links, and g is a coupling constant (identified with l/rn^ from eq. (3.1)) that has been tuned to 
a particular critical value that is 0(1). Unfortunately, little seems to be known about products of 
random matrices beyond dimension two [16]. Therefore we analyze instead a toy model where we 
define a "correlator" Q as a product of random numbers, and an "energy" = lim^^oo where: 



1 



Cz^YK^+gcpi) , = — ln(C 



(4.2) 



where < g < I and the (pi are independent and identically distributed random numbers with a 
uniform distribution on the interval [—1,1]. The exact value for the energy is obviously = 
since the statistical average of the correlator is (Q) = 1. The cumulants of the variable Y = In(CT) 
are given by 



ilog(l-g2) + 



tanh \g) 



1 



(-ir 



Li 



l-n 



(2tanh-^(g))^ 



V-gJ 



for n > 2; for g < 1 one finds that the Kn/n\ rapidly decrease as n increases. 

In Fig. 6 we show the results of a simulation where we compute S'^ for g -- 
1,...,1000. At each value of T we independently generated an ensemble of values for Q of 



and T 
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Figure 5: The quantities — ^ || and ^ as a function of N for unitary fermions at late times on a lattice 
of size L = 10, compared to mean field prediction (dashed lines). 
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Umax 2 



50 q- 250 500 1000 

Figure 6: Simulation of the energy (o^ for the toy model eq. (4.2) with g = ^. The exact answer is S'^ = 
(black line); exact values of the expansion eq. (4.1) truncated at order ^max = 2, 3 are indicated. 



Method 




Stat, error 


syst. error 


conventional 


0.014932 


0.002485 






-0.002159 


0.000304 


-0.002165 


Kn<3 


-0.000412 


0.001618 


-0.000324 


K«<4 


-0.000647 


0.008379 


0.000050 


K«<5 


-0.001794 


0.037561 


3.34 X 10"^ 


Km<6 


0.010943 


0.147739 


-1.22 X 10"^ 



Table 1: (o determined from 250 blocks of 50,000 configurations each for the model with t = 1000, 
8=1/2. 

size N = 50,000. From that ensemble we computed S'x by (i) using the conventional estimator 
= — ^ InCr (blue), which shows a striking systematic error for t > 50, and statistical noise in- 
creasing up to T 500 but decreasing beyond that; (ii) using eq. (4.1) truncated at rimax = 2 using 
sample cumulants Kn (green), showing a T-independent systematic error with smaller but slowly 
growing statistical error; (iii) eq. (4.1) truncated at rimax = 3 (red) with a negligible constant sys- 
tematic error but a larger statistical error. Evidently, one trades systematic error for statistical error 
by truncating eq. (4.1) at increasingly large Hmax- 

Table 1 displays results of a simulation of 1.25 x 10^ configurations blocked into 250 blocks 
of 50,000 each, for the model eq. (4.2) at T = 1000 and g = 1/2. We give the conventional 
estimate = — l/rlnC^ and estimates based on the cumulant expansion eq. (4.1) truncated at 
various rimax, where the exact value is ^ = 0. For each method we give the computed value for co 
and the statistical error; for the cumulant expansion we also give the exact systematic error from 
truncating eq. (4.1) at n = rimax using our analytic expressions for K^. These numbers show how 
the conventional method gives a wrong answer with deceptively small statistical error. For the 
cumulant expansion one sees again the trade of systematic error for statistical error with increasing 
rimax' Table 1 shows the combined error is minimized for rimax = 3, justified by noting that the 
^max = 4 result with statistical errors encompasses the rimax = 3 result; we suggest this as a practical 
algorithm for determining where to truncate the cumulant expansion in general. 

Leaving the toy model and returning to real simulations of unitary fermions. Fig. 7 shows how 
the cumulant expansion works for 50 trapped unitary fermions [11], where truncating the expansion 



0.01 
0. 
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Figure 7: Energy for 50 unitary fermions in a harmonic trap of frequency CO, 10^ configurations; fits per- 
formed using truncated at order nmax- Inset: conventional effective mass meff{t) = logC(T)/C(T+ 1) 
(green) and corresponding fitted cumulant effective mass with nmax = 4 (blue). 
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Figure 8: Effective mass plot for the same data displayed in Fig. 3, including ^ = 2,3,4,5 terms in the 
expansion eq. (4.1). The conventional analysis is given by gray triangles; the cumulant expansion results 
are blue circles. In each case we have marked the plateaus we extract, with bands whose widths signify the 
computed error. If the correlator distribution was exactly log-normal, the ^ = 2 analysis would be exact. 



at Umax = 4 is supported by the data. In Fig. 8 we show the improvement to the same data displayed 
above in Fig. 3 when we use the expansion eq. (4.1) to order n = 2,3,4,5. Again, as in the toy 
model, it is evident from this plot how increasing the number of cumulants included in the truncated 
version of eq. (4.1) decreases the systematic error, but at the cost of increased statistical error. 

5. Discussion 

Heavy-tail distributions are likely to be ubiquitous in A/^-body simulations, including multi- 
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Figure 9: Statistical distribution of the log of the real part of AA correlators at different temporal extent t 
suggesting a nearly log-normal distribution for the correlators. Each curve represents 100,000 samples; a 
small number of correlators with negative real part were discarded. Figure from W. Detmold (NPLQCD). 



baryon computations in lattice QCD. An example of the latter is given in Fig. 9, consisting of 
histograms of the log of the AA correlator, as computed by the NPLQCD collaboration from lattice 
QCD simulations. This plot shows evidence that the correlator must have a nearly log-normal 
distribution during the times considered. Presumably at very late times the distribution of the 
correlator will go over to a Gaussian with large variance and small mean, following the Lepage 
argument; however, at these earlier times that argument evidently fails, and the types of statistical 
techniques described here and in Ref. [1] might improve the analysis of such QCD data. It is 
interesting to speculate whether evidence for log-normal behavior in QCD implies that a mean 
field description could be valid for QCD correlators at moderate times, as discussed in §3; note the 
similarity between the QCD example Fig. 9, and the analogous plot for unitary fermions. Fig. 4 
— as well as the different way the two distributions scale with time. One might expect that if 
such a mean field description is approximately valid for multibaryon correlators in QCD, as baryon 
number is increased that description would become even more accurate, and log-normal behavior 
would become ubiquitous. 

It is also interesting to wonder whether the expansion in eq. (4.1) could be reformulated in an 
effective field theory language, where higher cumulants play the role of irrelevant operators. The 
reason the correlator distribution is flowing toward a simple universal form (log-normal) is presum- 
ably because fermions propagating through a single background field configuration are sampling 
many values of that field which are nearly statistically independent. What one measures is a sum of 
products of weakly correlated random matrices. The log of the correlator is in some sense averag- 
ing over short range fluctuations in the field, and it must be this averaging process — similar to RG 
blocking — which accounts for the flow of the distribution toward a Gaussian infrared attractor, 
with deviations characterized by the cumulants of the log of the correlator. Analogies between the 
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RG and statistics are well known — an example being the Central Limit Theorem (see, for example, 
[14]), or extreme statistics [15]. If the physical basis for attractors for the statistical distribution of 
correlators was better understood (not much is known about the distributions of products of random 
matrices, let alone sums of such products, but see [16]), including a systematic parameterization of 
deviations from the fixed point — analogous to the systematic inclusion of the effects of irrelevant 
operators in an effective field theory — then it is conceivable that one could greatly enhance one's 
ability to find a signal hiding in what superficially looks like useless noise. 
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